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SECTION  I 
INTRODUCTION 


In  recent  years,  significant  progress  has  been  made  toward  developing 
a  useful  method  for  predicting  steady  transonic  airloads  over  airfoils  and  to 
some  degree  for  finite  wings.  Despite  this  progress,  however,  very  few 
satisfactory  methods  exist  for  calculating  unsteady  transonic  flows,  as  evi¬ 
denced  in  a  recent  survey  paper  by  Bland  (Ref.  1).  Only  recently  have  num¬ 
erical  solutions  via  the  finite  difference  relaxation  techniques  been 
presented  for  two-dimensional  airfoils  executing  harmonic  motion.  As  an 
alternative  solution  technique,  the  authors  have  recently  investigated  the 
application  of  the  finite  element  technique  for  solution  of  the  transonic  flow 
problem,  and  computed  successfully  transonic  flow  over  nonlifting  thin  air¬ 
foils  (Ref.  2).  Compared  with  the  general  finite  difference  relaxation  tech¬ 
niques,  the  finite  element  method  provides  a  more  flexible  mesh  arrangement, 
element  shape  and  size  to  cope  more  effectively  with  complex  geometry  and 
boundary  conditions.  Also,  higher  order  approximations  can  be  readily  imple 
mented  to  increase  computational  efficiency. 

In  the  present  study,  the  finite  element  technique  is  extended  to  com¬ 
pute  lifting,  unsteady  transonic  flows  over  airfoils.  Since  the  objective  of 
the  present  study  is  to  develop  an  efficient  and  accurate  numerical  algorithm 
for  the  analysis  of  unsteady  transonic  flow  over  thin  airfoils,  the  present 
formulation  is  therefore  based  on  the  small  disturbance  but  nonlinear  tran¬ 
sonic  potential  equation  for  inviscid,  compressible  fluid.  Unlike  most  exist¬ 
ing  techniques  (Refs.  3  through  6)  for  solving  the  small  disturbance  potential 
equation,  the  present  approach  is  aimed  at  solving  directly  the  unsteady 
transonic  equation,  with  both  time  derivative  terms  retained.  Thus  the 
present  algorithm  can  be  applied  to  compute  a  much  wider  class  of  transonic 
flow  problems,  including  steady,  oscillatory  or  transient  solutions,  either 
with  or  without  angle  of  attack.  For  oscillatory  flow,  no  assumption  is  made 
regarding  the  oscillating  frequencies,  nor  must  the  unsteady  perturbation  be 
small  compared  to  the  mean  steady  solution. 
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On  the  numerical  aspects,  the  present  numerical  algorithm  is  developed 
using  the  cone  opt  of  finite  elements  in  conjunction  with  the  method  of  weighted 
residuals.  However,  to  treat  properly  the  mixed  flow  behavior,  the  residuals 
are  modified  by  introducing  a  "one-sided  assembly"  scheme  in  the  supersonic 
region.  Also,  a  patching  technique  has  been  developed  to  match  the  finite 
element  solution  constructed  in  a  moderately  large  domain  with  an  asymptotic 
solution  for  the  far  field  to  avoid  having  to  use  a  very  large  domain  otherwise 
needed  in  computations.  The  basic  element  presently  used  is  a  product  of  an 
element  in  space  and  an  element  in  time.  The  element  in  space  has  a  cubic 
expansion  inside  the  element,  with  nodal  unknowns  representing  the  per¬ 
turbed  velocity  potential  and  the  two  perturbed  velocity  components;  while  the 
element  in  time  is  a  quadratic  Lagrange  element.  The  resulting  system  of 
algebraic  equations  is  banded,  although  nonsymmetric,  which  is  solved  con¬ 
veniently  by  Gaussian  elimination. 

The  computations  of  steady  transonic  flow  over  lifting  airfoils  are  de¬ 
scribed  in  Section  II,  including  relevant  equations,  the  numerical  approach  and 
computed  results.  Considered  in  Section  III  are  the  calculations  of  unsteady 
transonic  flow,  which  indeed  can  be  applied  also  to  compute  steady  transonic 
flow.  The  final  section  contains  conclusions  and  recommendations  for  future 
studies . 


SECTION  II 

STEADY  TRANSONIC  FLOW  OVER  LIFTING  AIR  FOILS 


The  small  perturbation  theory  for  steady  transonic  flow  over  lifting 
airfoils  and  the  numerical  approach  for  solving  the  nonlinear  transonic 
equation  are  discussed  in  this  section.  First,  the  theory  including  the  small 
perturbation  potential  equation  with  associated  boundary  conditions,  together 
with  related  secondary  unknowns,  are  summarized.  A  finite  element  algo¬ 
rithm  based  on  a  modified  least  squares  method  of  weighted  residuals  is 
then  described  regarding  the  finite  element  formulation,  special  considera¬ 
tions  of  supersonic  regions,  the  imposition  of  boundary  conditions,  and  itera¬ 
tive  procedures  used  to  solve  the  resulting  system  of  nonlinear  algebraic 
equations.  Finally,  typical  flow  fields  for  a  6%  thick  circular  arc  and  a 
NACA  64A  410  airfoil  are  computed  to  demonstrate  the  feasibility  and  appli¬ 
cability  of  the  present  approach. 

1.  ASSUMPTIONS  AND  BASIC  EQUATIONS 

Since  the  objective  of  the  present  study  is  to  compute  transonic  flows 
over  thin  airfoils,  small  disturbance  theory  is  generally  adequate  and  hence 
used  as  the  basis  of  this  numerical  model.  With  this  theory,  perturbation  of 
the  flow  due  to  the  presence  of  airfoils  is  assumed  to  be  everywhere  small, 
the  embedded  shock  waves  are  assumed  to  be  weak  and  boundary  layer  effects 
are  neglected.  Also,  the  flow  is  assumed  to  be  isentropic  and  inviscid.  These 
assumptions  lead  to  a  nonlinear  equation  of  mixed  elliptic -hyperbolic  type  for 
general  two-dimensional  steady  transonic  flow  problems. 
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Boundary  Conditions 

v  -  (1  +  u)  3S-  s  0 
cix 

\4>  =  0  at  infinit  y  (3) 

where  =  perturbed  velocity  potential  function,  u,  v  =  perturbed  velocity 

components  in  the  x-  and  y -directions ,  respectively,  M  =  freestream  Mach 

00 

number,  y  =  ratio  of  specific  heats,  taken  to  be  1.4  for  air,  and  g  =  function 
of  x  defining  the  geometry  of  the  airfoil  including  angle  of  attack.  Equation  (11 
is  in  dimensionless  form  and  the  x-axis  is  aligned  with  the  undisturbed  flow 
direction.  The  dimensional  (with  primes)  and  nondimensional  quantities  are 
related  by 

x '  y*  ,  ,  _  0 ' 

x  "  T’  y  =  C  and  <t>  -tTT- 

00 

where  c  and  U  are  the  characteristic  length  and  speed,  which  are  currently 

00  a* 

taken  as  the  chord  length  and  the  freestream  speed,  respectively. 

For  flow  over  lifting  airfoils,  due  to  the  presence  of  circulation,  the 
velocity  potential  function  is  not  continuous  but  possesses  a  jump  equal  to  the 
circulation  strength.  For  this  reason,  a  branch  cut  must  be  made  to  allow 
the  potential  function  to  acquire  a  jump,  but  with  continuous  velocities. 
Therefore,  along  a  branch  cut,  the  following  conditions  must  be  imposed. 

*+  -  =  r 

u+  -  u  =  0  (4) 

v+  -  v  =0 

in  which  "+"  and  designate  the  upper  and  lower  surfaces  of  the  branch 
cut,  and  F  represents  the  strength  of  circulation. 


Once  the  flowfield  solution  in  terms  of  the  perturbed  velocity  potential 
has  been  obtained,  all  secondary  unknowns  can  then  be  calculated.  These 
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In  the  above,  U  =1,  the  normalized  freestream  speed,  a  =  local  speed  of 
sound,  M  =  loc^  Mach  number,  V  =  the  flow  speed,  p  =  local  static  pres- 
p^  =  stagnation  pressure,  and  Cp  =  pressure  coefficient. 


sure, 


2.  NUMERICAL  APPROACH 

The  finite  element  method,  in  conjunction  with  the  least  squares  method 
of  weighted  residuals,  is  used  herein  to  solve  numerically  the  small  pertur¬ 
bation  transonic  equation.  Summarized  in  this  subsection  are  the  finite  ele¬ 
ment  formulation,  considerations  of  the  embedded  supersonic  regions, 
imposition  of  boundary  conditions,  and  finally  iterative  procedures  used  to 
solve  the  resulting  system  of  nonlinear  algebraic  equations.  The  present 
numerical  approach  is  a  direct  extension  of  a  previous  study  for  transonic 
flow  over  nonlifting  airfoils.  More  detailed  discussion  is  given  in  Ref.  7. 
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I'  inito  Element  Formulation 


To  solve  the  present  problem  by  finite  element  methods  in  conjunction 
with  the  least  squares  approach,  a  set  of  locally  defined  trial  functions  with 
undetermined  parameters  is  first  assumed  as  the  approximate  solution,  and 
the  integral  expression  for  the  square  of  errors  committed  by  the  approxi¬ 
mate  solution  is  formulated.  Then  the  integral  of  square  errors  in  the  entire 
domain  is  minimized  with  respect  to  the  undetermined  parameters  to  yield  a 
system  of  algebraic  equations.  In  actual  computations,  the  minimization 
process  is  performed  at  element  level  and  then  an  assembling  process  is 
invoked  to  obtain  the  system  of  algebraic  equations. 


Written  in  the  usual  manner  with  repeated  indices  implying  summation, 
the  approximate  solution  has  the  following  form 

4  =  N.  4-  (9) 

V  1  T  l 


in  which  NF  s  represent  the  shape  functions  and  ^.'s  are  the  undetermined 
parameters  at  nodal  points.  The  residual  then  becomes 
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from  which  an  integral  expression  for  the  square  errors  is  obtained  as 
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--/ft dA 


(11) 


Upon  the  minimization  of  X  with  respect  to  the  undetermined  parameters,  a 
system  of  algebraic  equations  in  the  following  form  is  obtained 


S..  4  =0  (12) 
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Herein  the  banded  system  matrix  S.  .  is  defined  as 

U 


S..  =  ff P.  Q.  dA 

lJ  JJ  1  J 


with  P.  and  Q.  defined  by  the  following  expressions 
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As  stated  earlier,  the  system  matrix  is  obtained  by  combining  appropri¬ 
ate  contributions  from  all  the  elements.  The  element  matrices,  in  turn,  are 
evaluated  effectively  by  numerical  integration  to  avoid  the  tedious  and  error 
prone  algebraic  manipulations.  It  is  to  be  noted  that  although  the  system  of 
equations,  Eq .  (12),  is  homogeneous,  it  does  possess  a  non -trivial  solution 
once  the  boundary  conditions  are  imposed. 


The  basic  element  used  in  the  present  program  is  the  nonconforming 
cubic  triangular  element  developed  by  Bazeley  et  al.  (Ref.  8).  Also  used  in 
the  program  are  quadrilateral  and  trapezoidal  elements  constructed  from 
these  triangular  elements.  These  three  types  of  elements  can  be  mixed  and 
used  freely  in  the  entire  flow  region  except  that  only  trapezoids  should  be 
used  to  cover  adequately  the  supersonic  region  in  order  to  enact  the  special 
treatment  discussed  later. 


The  basic  triangular  element  is  shown  in  Fig.  1,  which  at  each  vertex 
has  the  function  itself  and  its  two  first  derivatives  (velocity  components)  as 
undetermined  parameters.  This  type  of  element  was  adapted  mainly  because 
boundary  conditions  of  both  Dirichlet  and  Neumann  types  can  be  imposed  with 
equal  convenience.  In  addition,  because  velocities  at  nodes  are  treated  as 
primary  unknowns,  secondary  unknowns  such  as  local  Mach  number,  and 
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pressure  coefficients,  etc.,  can  be  computed  directly  without  resorting  to 
numerical  differentiation,  thus  assuring  higher  accuracy.  Furthermore, 
the  use  of  higher  order  elements  can  usually  improve  computational  effi¬ 
ciency  as  evidenced  in  most  finite  element  analyses. 

As  is  known,  the  elements  presently  used  are  continuous  only  at 
nodal  points  but  not  across  element  boundaries,  which  apparently  violates  C* 
continuity  between  elements  as  normally  required.  However,  our  experience 
as  well  as  applications  by  others  do  not  seem  to  indicate  that  such  a  require¬ 
ment  is  necessary,  though  it  is  obviously  a  sufficient  condition  for  monotonic 
convergence  as  the  element  mesh  is  refined.  For  the  problems  of  plate  bending, 
it  was  concluded  that  the  condition  ecessary  to  guarantee  convergence  is 
that  the  elemert  used  should  be  capable  of  representing  constant  curvatures 
(second  derivatives)  within  the  element,  as  evidenced  by  several  cases  of 
plate  bending  analysis  (Ref.  8).  The  inter- element  continuity  requirement 
for  solving  the  present  second  order  equation  using  least  squares  should  be 
the  same  as  that  for  solving  a  fourth  order  equation  using  the  variational 
principle,  as  the  integral  forms  for  both  problems  involve  up  to  second 
derivatives  only.  In  the  present  study,  use  has  been  made  of  these  simpler 
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"non-conforming"  elements  and  results  appear  to  be  adequately  accurate. 

It  is  to  be  emphasized  that  "conforming"  elements  (with  continuous  normal 
derivatives  between  elements)  do  exist  and  can  be-  readily  implemented.  With 
these  elements,  however,  more  sophisticated  techniques  for  handling  the  dis¬ 
continuity  along  the  path  of  shock  wave  may  also  become  more  desirable.  Some 
studies  on  the  effects  of  "non-conforming"  versus  "c onforming "  el ements  ,  to¬ 
gether  with  more  sophisticated  schemes  for  handling  shocks,  are  indeed  highly 
desirable  and  useful. 

Considerations  of  Supersonic  Regions 

As  is  well  known,  the  finite  element  method  is  a  powerful  tool  for  solv¬ 
ing  problems  governed  by  elliptic  equations.  However,  for  problems  governed 
by  either  parabolic  or  hyperbolic  equations,  the  conventional  finite  element 
assembly  technique  must  be  modified  so  that  the  proper  influence  of  solution 
in  time  (or  time-like  direction)  is  considered.  For  the  present  problem,  the 
x-axis  is  the  time-like  direction  which  implies  that  the  solution  at  the  upwind 
station  will  determine  the  solution  at  its  downwind  station  but  not  the  contrary. 
This  consideration  leads  to  the  well  known  upwind  influence  finite  difference 
operator  (backward  finite  difference).  An  equivalent  technique  in  finite  element 
analysis  has  been  devised  and  applied  to  solve  the  transonic  flow  equation. 

Consider  the  rectangular  element  sketched  below  with  upwind  station  I 
and  downwind  station  II,  each  having  two  nodal  points.  With  the  element  type 
chosen,  the  element  matrices  can  be  constructed  in  the  usual  manner.  How¬ 
ever,  before  assembling  the  element  matrices  into  the  system  matrix,  the 
coefficient  of  the 
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first  term  in  the  transonic  equation  should  be  evaluated 

C  =  1  -  M2  -  M2  (1  +  y)  u 
oo  oo 


(16) 


The  sign  of  the  coefficient  C  being  positive,  zero,  or  negative  will  define  the 
equation  as  elliptic,  parabolic  or  hyperbolic.  If  C  is  non-positive  for  all 
nodes  in  the  element,  the  rows  representing  the  ;mproper  downwind  influence 
on  the  solution  at  an  upwind  station  are  ignored  during  assembly.  This  feature 
is  taken  care  of  automatically  in  the  program,  requiring  only  a  little  care  with 
the  nodal  arrangement  in  the  element.  In  the  anticipated  supersonic  region, 
element  node  points  should  be  arranged  in  the  order  as  indicated  in  the  above 
sketch,  starting  with  the  upper  left  corner  node  and  proceeding  in  a  counter¬ 
clockwise  direction.  On  the  other  hand,  if  the  sign  of  C  is  positive  at  any  of 
the  four  nodes,  no  special  assembling  is  invoked.  As  stated  earlier,  only 
trapezoidal  elements  are  to  be  used  in  the  anticipated  supersonic  and  mixed 
flow  region,  while  triangular  elements  and  quadrilaterals  consisting  of  only 
two  triangles  are  considered  unsuitable  because  of  their  bias  nature.  How¬ 
ever,  these  two  latter  types  of  elements  can  be  used  effectively  in  the  subsonic 
flow  region. 

Imposition  of  Boundary  Conditions 

As  stated  earlier,  the  imposition  of  boundary  conditions  for  the  present 
problem  can  be  carried  out  conveniently  because  the  elements  presently  used 
have  the  function  value  and  the  two  first  derivatives  as  primary  unknowns.  The 
associated  boundary  conditions  thus  can  be  treated  as  the  es sential  t  ype,  i.e., 
having  prescribed  values.  Standard  finite  element  methodology  is  therefore 
followed  by  assembling  first  an  unconstrained  problem  and  then  modifying  the 
matrix  equations  accordingly. 


On  the  airfoil,  the  nonlinear  form  of  flow  tangency  on  the  airfoil  is  im¬ 
posed  by  replacing  the  algebraic  equation  for  v.  by 
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dx 
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in  which  and  are  unknown  parameters  at  node  "i"  on  the  airfoil.  When 
the  linear  form  is  desired,  the  corresponding  equation  becomes 


v. 

i 


= 

dx 


(18) 


and  is  applied  along  the  chordline. 

For  the  far  field  boundary  condition,  Eq.  (3)  is  valid  only  at  infinity, 
and  thus  an  infinite  domain  is  required  in  computations  if  the  condition  is  to 
be  imposed  directly.  Although  an  approximation  can  be  made  by  assuming  Eq.  (3) 
to  be  valid  on  a  boundary  at  a  finite  distance  from  the  airfoil,  a  large  computa¬ 
tional  network  is  usually  required  to  obtain  a  solution  with  adequate  accuracy 
due  to  the  fact  that  for  lifting  airfoils  at  transonic  speeds,  the  disturbances 
created  by  the  body  decay  very  slowly.  There  are  several  ways  to  cope 
with  this  problem  so  that  computations  can  be  restricted  to  a  finite  region. 

One  approach  is  to  map  the  region  exterior  to  the  airfoil  into  the  interior  of 
a  circle  upon  which  the  computations  are  performed  and  the  final  solution  is  ob¬ 
tained  by  a  subsequent  inversion  (Ref,  9).  This  technique  is  almost  ideal  for 
computations  involving  a  single  airoil,  but  its  extension  to  three-dimensional 
cases  or  flow  over  multiple  airfoils  is  not  evident.  Another  approach  is  to 
match  the  numerical  solution  of  the  near  field  with  an  analytic  representation 
for  the  far  field  to  satisfy,  in  effect,  the  far  field  boundary  conditions,  as  em¬ 
ployed  in  finite  difference  relaxation  techniques  (Refs.  10  and  11).  In  doing  so, 
computations  can  be  performed  in  the  physical  plane  by  using  a  finite  domain 
without  sacrificing  too  much  accuracy  and  more  importantly,  it  provides  a 
straightforward  extension  to  treat  three-dimensional  and  multiple  airfoil 
problems.  This  concept  is  pursued  with  appropriate  modifications  in  the 
present  finite  element  solution  technique  so  that  the  circulation  and  hence 
lifting  force  is  computed  by  a  systematic  approach. 

Figure  2  depicts  a  lifting  airfoil  under  free-flight  conditions.  In  the 
analysis,  a  moderately  large  domain,  R',  is  taken  to  construct  the  finite  ele¬ 
ment  solution,  while  an  asymptotic  solution  with  undetermined  parameters 


11 


R":  Far-Fieid  Solution 


Finite  Element  Solution 


S,:  Matching  of  Solutions 


Matching  of  Localized  Finite  Element  Solution 
with  Far-Field  Solution 


is  assumed  to  be  valid  in  the  far  field,  R."  The  two  solutions  are  then  matched 
along  the  common  boundary,  Sj,  along  which  the  unknowns  are  then  expressed 
in  terms  of  the  far-field  undetermined  parameters.  Along  the  branch  cut,  S^.  th 
conditions  of  the  velocity  being  continuous  but  the  velocity  potential  possessing 
a  jump  equal  to  the  circulation,  r>  are  imposed  iteratively.  Through  these 
procedures,  the  far-field  parameters  will  be  determined  systematically  rather 
than  merely  by  trial  and  error. 


The  far-field  potentials  for  both  lifting  airfoils  and  three-dimensional 
wings  at  transonic  speeds  were  developed  by  Klunker  (Ref.  12).  The  expression 
for  a  two-dimensional  airfoil  is  currently  used  and  has  the  following  form 
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The  first  term  on  the  right-hand  side  corresponds  to  a  free  vortex  and  repre¬ 
sents  the  lifting  effects;  the  second  term  corresponds  to  a  source  distribution 
whose  strength  is  related  to  the  airfoil  geometry  g(x);  and  the  third  term, 


which  arises  from  the  nonlinearity  of  the  flow  equations,  has  the  form  of  a 

2 


doublet  with  its  strength  given  by  the  local  value  of  u  and  is  to  be  integrated 
over  the  flow  domain  under  consideration.  Klunker  also  showed  that  the 
contributions  from  the  thickness  and  the  area  integrals  are  of  higher  order 
effects.  For  these  reasons,  only  the  first  term  representing  the  lifting  effects 


is  used  currently  for  the  far  field.  The  solution  along  the  outer  boundary  is 


updated  systematically,  using  the  computed  potential  jump  at  the  trailing  edge. 
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The  boundary  condition,  Eq.  (4),  along  the  branch  cut  can  be  considered 
as  constraints  to  the  unknown  parameters.  These  constraints  can  be  imposed 
through  introducing  Lagrangian  multipliers.  As  an  alternative,  a  scheme 
without  the  need  of  introducing  Lagrangian  multipliers  was  devised  and  im¬ 
plemented.  The  scheme  is  essentially  based  on  the  concept  of  using  chain 
rule  in  minimizing  the  integral  of  square  errors.  More  specifically,  the 
algebraic  equatians  are  first  generated  assuming  no  constraint  on  the  original 
parameters.  However,  because  of  constraint  on  certain  unknown  parameters, 
the  algebraic  equations  corresponding  to  constrained  parameters  are  then 
manipulated  according  to  chain  rule  and  the  specific  constraint  to  yield  one 
equation  less  for  each  constraint.  And  finally,  one  unknown  parameter  could 
be  eliminated  for  each  constraint,  thus  making  the  total  number  of  equations 
and  unknowns  equal.  A  variation  of  this  scheme  is  to  insert,  in  place  of  the 
equation  eliminated,  the  constraint  and  leave  the  total  number  of  unknown 
parameters  unchanged.  The  latter  approach  has  the  advantage  of  simpler 
bookkeeping  and  was  therefore  implemented  in  the  computer  program.  For 
example,  to  impose  the  condition  u  =  u  for  a  pair  of  nodes  along  the  branch 
cut,  the  equation  originally  generated  for  u  is  first  added  to  the  equation  for 
u+,  and  then  in  its  place  the  constraint  equation  u+  -  u  =  0  is  inserted. 

Iterative  Procedures 

With  the  equations  properly  assembled  and  boundary  conditions  imposed, 
the  system  of  nonlinear  algebraic  equations  is  finally  solved  by  iterative  pro¬ 
cedures  of  the  form 

S.  .  (f)  *jn)  =  L.  (20) 

to  solve  for  the  solution  i>  ^  at  the  nfc^  iteration.  The  function  i>  is,  in  turn, 
defined  as 

+  =.  9  +  (1  -  9)  (21) 
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in  which  9  is  a  relaxation  factor  in  the  range  0  <  9  <1.  For  subsonic  flow, 

9  is  usually  taken  as  unity;  for  flow  in  the  transonic  regime,  however,  it  has 
been  found  that  the  use  of  an  under-relaxation  factor  is  required.  The  value  of 
9  should  decrease  from  unity  for  barely  critical  flow  to  0.4  or  even  0.3  for 
supercritical  flow.  The  optimum  under  -  relaxation  factor  generally  depends 
on  the  freestream  Mach  number  and  the  mesh  being  used.  At  the  present,  no 
systematic  investigations  have  been  done  in  this  regard,  and  therefore  its  value 
should  be  selected  by  numerical  experimentations.  Generally  speaking,  a 
smaller  under -relaxation  factor  will  make  the  solution  more  stable  but  at  the 
same  time  tends  to  slow  the  rate  of  convergence. 


Equation  (20)  is  to  be  solved  subject  to  certain  prescribed  convergence 
criteria.  The  one  presently  used  is  that  the  relative  change  of  local  Mach 
number  between  two  consecutive  iterations  should  be  less  than  a  prescribed 
small  number  for  all  the  nodes  in  the  flow  field,  that  is 


M<">  - 
M(n> 


<  e 


(22) 


The  value  of  e  in  most  umerical  computations  is  in  the  range  0.001  <  e  < 
0.005. 


The  numerical  procedures  can  thus  be  summarized  as  the  following; 


1.  Construct  finite  element  algebraic  equations  for  the  inner  field 
based  on  the  method  of  weighted  residuals,  presently  with  the 
least  squares  approach  with  modifications  made  in  supersonic  regions. 


2.  Apply  boundary  condition  on  the  airfoil  to  satisfy  flow  tangency 
condition  at  nodal  points. 


3.  Impose  the  continuity  of  velocity  components  along  the  branch 
cut  as  described  earlier. 


4.  Specify  the  parameters  along  the  common  boundary  of  inner  and 
outer  field  in  terms  of  the  value  of  circulation  strength. 


15 


r>.  Solve  for  the  flow  field  and  hence  the  jump  in  potential  at  the 
trailing  edge,  which  is  then  used  as  an  improved  solution  in 
updating  the  finite  element  algebraic  equations  and  the  circula¬ 
tion  along  the  outer  boundary,  respectively. 

6.  Repeat  tne  above  procedures  until  a  convergent  solution  is  ob¬ 
tained,  such  as  until  the  relative  change  of  Mach  number  falls 
within  a  small  range  of  tolerance  or  the  jump  in  potential  function 
across  the  branch  cut  is  essentially  a  constant. 

3.  NUMERICAL  RESULTS 

To  demonstrate  the  applicability  of  the  present  approach,  flow  fields 
over  a  6%  thick  circular  arc  and  a  NACA  64  A410  airfoil  have  been  computed. 
Results  obtained  by  the  present  approach  appear  to  compare,  in  general,  very 
well  with  experimental  data  and  those  obtained  by  relaxation  techniques,  while 
each  case  calculated  requires  only  a  few  minutes  of  CPU  time  on  a  Univac 
1108  computer . 

For  the  6%  thick  circular  arc,  calculations  were  made  for  cases  with 
angle  of  attack  equal  to  1  and  2  deg.  The  mesh  presently  used  is  shown 
in  Fig.  3,  with  202  nodes.  As  is  seen  in  the  figure,  the  mesh  arrange¬ 
ment  is  completely  flexible,  except  in  the  expected  supersonic  region  where 
i  rather  regular  mesh  is  needed  in  order  to  consider  properly  the  problem  of 
region  of  dependence.  The  predicted  results  for  the  two  cases  are  shown  in 
Mgs.  4  and  5,  respectively,  where  comparisons  are  made  with  experimental 
data  obtained  by  Knechtel  (Ref.  13).  It  is  seen  that  the  agreement  between  the 
predicted  results  and  experimental  data  is  excellent  except  near  the  leading 
and  trailing  edges.  These  noticeable  discrepancies  are  believed  mainly  due 
to  the  invalidity  of  the  small  perturbation  equation  in  these  regions.  In 
addition,  the  present  computations  are  based  on  inviscid  theory,  while  experi¬ 
mental  results  necessarily  involve  viscous  effect,  which  were  evidenced  in  the 
form  of  negative  aerodynamic  loadings  over  the  rear  portion  of  ths  airfoil. 

For  these  reasons,  the  predicted  results  are  considered  to  be  sufficiently 
accurate.  Among  all  of  the  cases  calculated,  the  greatest  difference  appears 
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to  be  near  the  leading  edge  for  the  critical  case  (a  =  2  ,  M^  =  0.855),  where 
a  large  peak  pressure  is  present  in  the  experiment.  However,  as  other 
cases  of  lower  and  higher  Mach  numbers  do  not  indicate  such  large  peak 
pressure  near  the  leading  edge  (see  Fig.7e  of  Ref.  13),  the  peaky  behavior 
for  the  present  case  is  somewhat  peculiar  and  questionable. 

To  test  the  present  approach  further,  the  flow  field  over  a  NACA  64 
A4  10  airfoil  was  also  computed.  This  airfoil  was  chosen  as  another  testing 
case  mainly  for  two  reasons:  (1)  the  airfoil  is  thicker  (with  10%  thickness 
ratio)  and  blunt-nosed,  thus  it  represents  a  more  severe  test  for  the  present 
method;  and  (2)  the  airfoil  has  been  studied  rather  extensively  and  both  ex¬ 
perimental  and  numerical  data  are  available  for  comparisons  (see  Ref.  11). 

To  date,  only  limited  computations  have  been  performed.  Preliminary  re¬ 
sults  are  depicted  in  Figs.  6,  7  and  8  for  cases  involving  various  Mach  numbers, 
either  with  or  without  angle  of  attack.  As  the  present  airfoil  has  a  blunt 
leading  edge,  the  boundary  condition  in  that  region  cannot  be  imposed  in  the 
usual  manner  by  specifying  the  vertical  velocity  component  as  v  =  dg/dx  =  oo . 
Several  numerical  experiments  have  been  conducted  to  cope  with  these  diffi¬ 
culties.  One  of  the  attempts  is  to  specify  the  normal  derivative  of  the  per¬ 
turbed  velocity  potential.  The  value  of  the  normal  derivative,  in  turn,  can 
be  determined  explicitly  from  the  flow  tangency  condition.  To  accomplish 
this  in  the  numerical  computations ,  the  normal  and  tangential  derivatives 
of  the  potential  function  are  used  as  unknowns  at  the  leading  edge,  instead  of 
the  derivatives  originally  in  the  x-  and  y-dir ections .  Thus  the  normal  de¬ 
rivative  of  the  perturbation  potential  can  be  conveniently  specified.  These 
thoughts  have  been  implemented  into  the  program.  However,  our  numerical 
experimentation  indicates  that  the  imposition  of  a  full  amount  of  perturbed 
normal  velocity  as  computed  from  flow  tangency  condition  (without  invoking 
the  small  perturbation  assumption)  does  not  produce  satisfactory  results  near 
the  leading  edge,  with  the  flow  in  that  region  being  impeded  excessibly.  The 
inaccuracy  of  the  solution  in  this  region  is  believed  partly  due  to  the  mesh  being 
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Experimental 
Finite  Difference 


O  Finite  Element  (I) 
A  Finite  Element  (II) 


b.  M  =  0.74,  a  =  0  deg 

(30 

6  -  Comparison  of  Chordwise  Pressure  Distribution  for  NACA 
64  A410  Airfoil  (a  =  0  deg) 


too  coarse  and  partly  due  to  the  fact  that  the  boundary  condition  used  does 
not  match  the  small  disturbance  equation  there.  In  an  attempt  to  make  the 
boundary  condition  at  the  leading  edge  compatible  with  the  small  perturbation 
equation  which  is  assumed  to  be  valid  there,  numerical  experimentations  were 
performed  by  specifying  only  a  fraction  of  the  perturbed  velocity  in  the  normal 
direction,  with  results  depicted  in  Figs.  6  and  7.  Method  I  corresponds  to 
one-fourth  of  the  normal  velocity  otherwise  required,  while  Method  II  has  a 
zero  normal  perturbed  velocity  specified.  As  is  seen,  the  pressure  distri¬ 
bution  near  the  leading  edge  is  quite  sensitive  to  the  boundary  condition  im¬ 
posed  there;  however,  the  comparison  for  the  remaining  portion  of  the  airfoil 
is  very  good.  Another  attempt  is  to  use  finite  values  of  slope  for  the  nodal 
points  at  the  leading  edge  and  impose  the  flow  tangency  condition  as  required 
by  Eq.  (17).  Figure  8  presents  results  obtained  by  using  at  the  leading  edge 
half  the  value  of  the  airfoil  slopes  at  0.5%.  Compared  with  the  predicted  re¬ 
sults  shown  in  Fig.  7,  some  improvement  near  the  leading  edge  is  apparent, 
but  the  expansion  on  the  upper  surface  is  still  not  properly  simulated.  Never¬ 
theless,  as  seen  in  the  figures  the  inaccuracy  near  the  leading  edge  is  local¬ 
ized  and  is  always  confined  to  the  first  few  nodal  points.  The  solution  accuracy 
could  be  improved  within  the  small  disturbance  assumption  to  a  certain  extent 
by  refining  the  mesh  near  the  leading  edge.  However,  because  the  core  storage 
capacity  normally  available  on  the  local  Univac  1108  computer  is  only  64K, 
some  additional  programming  effort  is  required  in  order  to  compute  the  prob¬ 
lem  with  a  much  finer  mesh;  therefore,  only  results  with  the  present  mesh  are 
presented. 
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SECTION  III 

UNSTEADY  TRANSONIC  FLOW 


In  this  section,  finite  element  procedures  are  described  for  solving 
the  problems  of  unsteady  transonic  flow.  Again,  the  formulation  is  based  on 
small  perturbation  theory.  However,  unlike  the  approach  developed  earlier 
for  NASA-Langley  (Ref.  6),  the  assumptions  of  harmonic  motion  and  the  un 
steady  disturbance  being  small  compared  to  the  mean  steady  solution  are 
removed  in  the  present  approach.  The  removal  of  these  assumptions  allows 
the  flow  field  to  be  computed  in  a  wider  range.  These  include  the  problems 
involving  transient  solution,  moderate  unsteadiness  and,  more  importantly, 
movement  of  shock  wave  location.  The  equations,  upon  which  the  present 
numerical  procedures  are  constructed,  are  summarized  in  Subsection  1.  Two 
numerical  approaches,  namely,  the  Galerkin  type  and  least  squares  method, 
which  have  been  investigated  to  integrate  directly  the  unsteady  transonic  equa¬ 
tion,  are  described  in  Subsection  2.  As  a  result  of  the  present  investigation, 
the  Galerkin  method  of  weighted  residuals  was  found  to  be  invalid  for 
computing  supercritical  flow,  thus  the  least  squares  approach  has  been  pur¬ 
sued.  To  date,  only  limited  computations  have  been  conducted  using  the  least 
squares  approach;  however,  results  obtained  thus  far  indicate  the  present  ap¬ 
proach  is,  in  general,  very  satisfactory.  Typical  results  are  summarized 
and  discussed  in  Subsection  3. 


1.  BASIC  EQUATIONS 

Based  on  small  perturbation  theory,  the  unsteady  transonic  flow  prob¬ 
lem  can  be  stated  as 


L<*>  = 


=  0 


(23) 


in  which 


a  =  a  +  b  d> 

,  x 
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with 


a  =  1  -  M 


oo 


b  =  -  M  (1  +  y) 
oo 


The  boundary  conditions  associated  with  Eq.  (23)  include: 


<1  Vanishing  of  Disturbance  at  the  Far  Field, 


V*  =  0 


(24) 


Flow  Tangency  Condition  on  the  Airfoil  Surface, 


DB 

Dt 


=  B  +  ( 1  +  ^  )  B  +  A  B  =0 

,t  r,x'  ,x  r ,  y  ,  y 


with  B  describing  the  instantaneous  airfoil  position  as 

B(x,y,t)  =  y  -  g(x)  -  6h(x,t)  =  0 

In  the  above,  g(x)  =  geometry  of  the  airfoil  at  mean  steady  position,  6  =  ampli¬ 
tude  of  oscillation,  and  h(x,t)  =  function  describing  the  airfoil  oscillation. 

Upon  substitution,  a  nonlinear  boundary  condition  is  obtained  as 


*.v  ‘  5<\x+V+e,x<1  +  *,*> 


(25) 


Which  is  to  be  imposed  on  the  mean  position  of  the  oscillating  airfoil,  i.e., 

B(x,  y,  t  )  =  0. 

7  mean 


•  Unsteady  Kutta  Condition  in  the  Wake 


Along  the  wake,  pressures  and  downwashes  are  equal  on  both  sides  of 
the  vortex  sheet.  Since  the  pressure  coefficient  is  given,  to  the  lowest  order, 
as 


2*  t  -  2 

I  1  >  -*■  9 
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the  linearised  boundary  conditions  in  the  wake  thus  become 


and 


+  ♦'* 


(26) 


(27) 


which  are  to  be  applied  along  the  upper  and  lower  surfaces  of  the  mean  wake 
position.  The  potential  values  on  the  upper  and  lower  surfaces  of  the  wake 
are  related  by 

+  T  (28) 


in  which  F  =  F(x,t)  for  unsteady  flows.  With  T  as  the  primary  unknown,  the 
condition  expressed  by  Eq.  (26)  can  also  be  written  as 


r  =  -  r  (26a) 

,t  ,x 

As  discussed  earlier,  for  lifting  airfoil  computations,  the  far  field  con¬ 
ditions  must  be  considered  properly.  To  do  this  an  asymptotic  expression 
with  undetermined  parameters  must  be  established  and  matched  with  the 
near -field  finite  element  solution.  In  the  present  study,  only  the  asymptotic 
solution  for  mean  steady  flow  has  been  established  and  implemented,  im¬ 
plying  that  the  unsteady  perturbation  at  far  field  is  small  compared  to  the 
mean  steady  solution. 


2.  NUMERICAL  APPROACHES 

The  finite  element  technique,  in  conjunction  with  the  Method  of  Weighted 
Residuals  (MWR),  is  used  to  solve  Eq.  (23)  with  its  associated  boundary  condi¬ 
tions.  The  approaches  undertaken  include  the  Galerkin  type  and  the  least 
squares  method.  At  this  stage,  results  indicate  that  the  Galerkin  approach 
can  be  used  to  compute  flow  in  the  subcritical  regime  but  does  not  seem  to 
be  valid  for  calculating  critical  transonic  flows.  The  least  squares  approach, 
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on  the  other  hand,  appears  to  be  satisfactory  when  critical  flow  is  considered. 
These  two  approaches  are  described  in  the  following. 

The  Galerkin  Approach 

The  Galerkin  approach  was  investigated  during  the  earlier  stage  of  this 
study  as  it  can  offer  certain  advantages  over  the  least  squares  formulation, 
especially  when  the  unsteady  transonic  equation  is  considered.  These  include; 

1.  The  formulation  by  Galerkin  criterion  is  noticeably  simpler  and 
will  require  less  computation  time. 

2.  With  appropriate  integration  by  parts  and  arrangement,  the  con¬ 
tinuity  requirement  between  elements  can  be  lowered  and  boundary 
conditions  can  be  treated  more  conveniently. 

With  the  Galerkin  approach,  an  approximate  solution  for  the  velocity 
potential  is  first  assumed  in  the  form 

<f>  -  N.(x,  y)$.(t)  (i  =  1  to  n)  (29) 

in  which  IT  =  shape  functions,  <j> ^  =  the  corresponding  undetermined  parameters, 
and  n  =  total  number  of  unknown  parameters.  Equation  (29)  is  then  substituted 
into  (23)  and,  with  hT  as  weighting  functions,  the  procedure  of  weighted  residuals 
is  applied  in  the  spatial  directions  to  obtain  a  system  of  Ordinary  Differential 
Equations  (ODE)  in  the  form 

A. .  4>\  +  B. .  4>.  +  C. .  <f>.  =  0  (30) 

J  iJ  J  lJ  J 

in  which  the  coefficient  matrices  are  banded  and  are  functions  of  space  only. 

To  integrate  Eq.  (30),  a  Galerkin  type  time  marching  scheme  is  used  to 
obtain  a  recurrence  relationship  which  is  implicit  in  nature.  For  instance, 
with  a  quadratic  expression  for  the  time  history  of  each  nodal  unknown  one  has 

=  Mk(t)  (k  =  1  to  3)  (31) 
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in  which  represents  the  shape  functions  in  time.  Again,  direct  substitu¬ 
tion  and  applying  the  method  of  weighted  residuals  (choosing  as  the 
weighting  function)  with  integration  carried  out  in  time  will  result  in  a  re¬ 
currence  equation  in  the  form 

(a.  A..  +  b.  B..  +  c.  C.  .)  ^  =  0  (32) 

k  ij  k  ij  k  ij  j 

where  the  vectors  a^,  b^,  c^  are  functions  of  time  only.  Equation  (32)  can  be 
rewritten  as 


(a,A..+b,  B..  +  c,  C..)^?  =-(a  A.  .  +  b  A.  .  +  c  C..)^m 
3  ij  3  ij  3  ij'  Tj  m  ij  m  ij  m  ij  J 


(33) 


(m  =  1  to  2) 


From  Eq.  (33),  the  solution  at  time  level  k  =  3  can  be  computed  from  solutions 
at  the  previous  time  steps.  This  process  will  continue  until  the  solution  sought 
is  obtained. 


A  program  based  on  the  above  formulation  has  been  developed  and  used 
first  to  solve  a  steady  problem  which  is  obtained  as  the  limiting  flow  for  large 
times  with  suitable  initial  conditions.  Presently  the  initial  conditions  are  ob¬ 
tained  by  placing  a  leaky  profile  in  the  desired  uniform  stream  and  then  im¬ 
pulsively  turn  off  the  leakiness  at  zero  time.  That  is,  the  flow  is  assumed  to 
be  uniform  for  t  <  0,  and  is  required  to  satisfy  the  tangential  boundary  condi¬ 
tion  on  the  airfoil  for  t  >  0.  With  a  6%  thick  circular  arc  airfoil,  convergent 

results  are  obtained  for  M  =  0.806  (subcritical)  and  M  =  0.84  (barely  crit- 

oo  oo 

ical).  However,  as  M  increases  to  0.86  (critical),  the  sequence  of  unsteady 

oo 

flows  oscillates  from  one  time  step  to  another  and  diverges.  This  phenomenon 
persists  whether  the  improper  downwind  influence  is  removed  or  not. 

In  attempting  to  resolve  the  above  problem,  several  numerical  experi¬ 
ments  (variations  of  the  Galerkin  approach)  were  then  carried  out  for  the 
steady  transonic  equation. 
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As  found  previously  (Ref.  2)  the  transonic  equation  should  not  be  cast 
into  an  equation  of  the  Poisson  type  and  treated  as  an  elliptic  problem.  Any 
numerical  scheme  so  constructed  would  inhere  in  it  the  improper  down¬ 
wind  influence  upon  upwind  solution  in  the  supersonic  region.  To  eliminate 
this,  it  is  more  convenient  to  apply  the  Oulerkin  approach  directly  without 
integration  by  parts  for  the  term  involving  ^  .  However,  in  order  to  obtain 

a.  better  conditioned  coefficient  matrix,  integration  by  parts  for  the  <(>  term 

>  y  y 

is  carried  out.  This  approach  was  also  found  unsuccessful  in  obtaining  a  con¬ 
vergent  solution  for  the  case  M  =  0.861. 

00 

Another  attempt  was  to  add  a  term  simulating  numerical  viscosity  for 
elements  in  the  supersonic  region.  The  term  added  is 

aAX$  (34) 

,  XXX 

In  the  finite  difference  approach,  AX  is  the  spacing  between  two  grid  points. 
When  the  transonic  equation  is  discretized  by  a  second  order  central  differ¬ 
ence  formula,  this  term  is  capable  of  eliminating  the  downwind  influence 
(Ref.  14).  In  our  formulation,  AX  is  taken  to  be  the  X  dimension  of  the  ele¬ 
ment.  Again,  no  satisfactory  results  were  obtained  for  the  case  of  M  =  0.861. 

00 

Other  possibilities  for  modifying  the  conventional  Galerkin  approach  to 
solve  the  transonic  equation  remain  to  be  exploited.  However,  work  in  that 
direction  could  be  time  consuming  and  expensive.  Based  on  its  success  with 
solving  the  steady  transonic  problem,  the  least  squares  approach  was  there¬ 
fore  investigated  to  solve  also  the  unsteady  transonic  equation. 

The  Least  Squares  Approach 

An  approximate  solution  for  the  perturbed  velocity  potential  can  be 
assumed  in  the  form 

}  =  N.  M.  <£k  (35) 

i  k  Ti 
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in  which  N  =  shape  functions  in  space,  M,  =  shape  functions  in  time,  and 

k  1  * 

4> ^  =  unknown  parameter  corresponding  to  node  i  and  time  level  k.  Shape  func¬ 
tions  in  space  generally  depend  on  the  particular  element  being  used  as 
discussed  in  Section  2.  The  shape  functions  in  time  presently  used  are  the 
quadratic  approximating  functions  shown  in  Fig.  9. 

Upon  substituting  Eq.  (35)  into  Eq.  (33),  the  residual  R  can  be  expressed 
in  terms  of  the  unknown  parameters  directly.  That  is, 

R  =  M.  (aN.  +  N .  )  -  (2  ML  tN.  +  M,  (36) 

l  k'  j,xx  j,  yy  00  k,t  J.X  k,tt 


where  a  is  approximated  by 


a  =  a  +  bM„  N  4> 
l  p,x  rp 


The  residual  can  be  rewritten  as 


where 


R  =  (A^  +  B^) 

J  J  J 


A .  =  M.  (a  N .  +  N.  ) 

J  k  j,  xx  j,  yy7 


Bk  =  -  M2  (2  M.  N .  +  M.  ..  N.) 

J  00  '  k.t  J,x  k,  tt  J7 


The  integral  expression  of 


Iff 


R  dx  dy  dt 


is  minimized  with  respect  to  the  unknown  parameters  at  time  level  k  =  3, 
with  intjgration  being  taken  over  the  physical  space  and  to  the  time  level 
k  =  3. 
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That  is 


3*; 


Ilf 


R  dx  dy  dt  =  0 


or 


Iff  i^Rd*iydt  =  ° 


(38) 


Letting 


W.  =  —  =  A3  +  B3  +  (M.  N.  <£k)  bM-  N. 
l  3  l  l  k  j,  xx  r  j  3  l,  x 


the  resulting  system  of  algebraic  equations  becomes 


///',! 


(39) 


M,  («N.  +  N .  )  -  M.  ( 2  M.  N .  +  M.  N .) 

k  j,xx  j,  yy'  kl  k,t  ],x  k,tt  j' 


dx  dy  dt  =  0 


(40) 


Equation  (40)  provides  us  with  a  recurrence  relation  to  solve  for  4>  ■  in  terms  of 

1  2  J 

$  .  and  .  . 

J  J 


A  program  based  on  Eq.  (40)  was  developed  and  used  to  check  out  first 
the  steady  state  solutions.  Through  these  numerical  experimentations,  it  was 
found  that,  in  order  to  stabilize  the  time  marching  scheme,  emphasis  must  be 
placed  on  the  last  time  step  for  terms  involving  space  derivatives  only,  i.e., 


R’  =  <“'Nj,,X+Nj.y>r»*j-<'2Mk,tNj,X+Mk,t«Ni»^  <4‘> 


with 


a'  =  a  +  b  N  $ 
P. x  P 


(42) 


instead  of  Eqs.  (36)  and  (37).  Such  modifications  are  motivated  by  observing 
the  implicit  finite  difference  scheme  for  the  wave  equation, 


tt 


xx 


=  0 


(43) 
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Take  a  pivoting  point  (i,  j),  (i  denotes  x  and  j  denotes  t),  assume  mesh  size 
in  x  is  k  and  in  t  is  h,  then  the  difference  equation  is 


+  *. 


izi 


"  2  ; ’+ 


-  \ 


2  'i  +  l 


i-  1 , 


1  - 


=  0 


(44) 


Through  Von  Newmann's  stability  analysis  by  Fourier  series  (Ref.  15),  this 

equation  can  be  shown  to  be  unconditionally  stable  for  all  h  and  k.  An 

immediate  observation  is  that  4  is  differenced  in  the  x  direction  at  the 

xx 

last  time  step.  Since  the  left-hand  side  of  the  finite  difference  equation  (44) 
is  in  fact  the  residual  at  a  discrete  node,  the  analogous  residual  for  the  pre¬ 
sent  transonic  equation  via  finite  element  approach  therefore  is  Eq.  (41). 

Note  that  M,  4  ^  in  Eq.  (36)  is  replaced  by  (M .  +  M_+M_)^>?in  Eq.  (41), 
k  J  3  l  c  3  j 

which  reduces  to  4 ■ ,  since  M,  +  M,  +  M.  =  1. 

J  12  3 

With  the  residual  established,  the  corresponding  weighting  functions 
can  be  readily  determined,  thus  leading  to  the  required  system  of  algebraic 
equations. 

It  remains  to  argue  that 

AR  =  R  -  R'  =  |M_  ( 4 ?  -  4?)  +  M.  ( 4 ?  -  N. 

[  2  'rj  V  1  V|  J,  yy 

+  |mi  (a*J  -  a’4j)  +  M2  (a^2  -  a'4p  (45) 

+  M,  (a  -  a')  4?  |n. 

3  J  J  J.  xx 

If  steady  flow  is  approached, 

43  =  4z  =  a  =  a'-=^  AR  =0  (46) 
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Therefore  this  approach  is  appropriate  for  computing  steady  flow.  In  this 
case,  the  role  of  AR  seems  to  be  the  same  as  the  numerical  viscosity  term 
introduced  in  the  implicit  finite  difference  scheme  through  truncation. 

Other  aspects  of  the  numerical  approach  are  similar  to  what  was  de¬ 
scribed  in  Section  2.2  for  steady  flow,  except  that  the  boundary  conditions 
imposed  along  the  branch  cut  are  noticeably  different.  Also,  a  process  of 
marching  in  time  is  taken,  instead  of  successive  iterations  used  on  the 
steady  counterpart. 

Along  the  branch  cut  (vortex  sheet  in  unsteady  flow),  Eqs.  (26)  and  (27) 
must  be  imposed  to  ensure  the  continuity  of  pressure  and  downwash.  To  im¬ 
pose  the  first  condition,  Eq.  (26),  for  a  pair  of  nodes  along  the  branch  cut,  the 
two  finite  element  equations  generated  originally  for  and  $  are  properly 

combined  to  yield  one  equation,  while  the  equation  for  <f)  +  is  replaced  by  the 
following  equation,  evaluated  at  time  level  t  =  t^: 

+(3)  (3)  (3)  (3)  (2)  (2)  (1)  (1) 

*,x  -  +M3,t(^  )=M2>t(*'  -4  )+Ml  t(*'  -4>+  )  (47) 

Where  M  is  the  shape  function  in  time  and  the  superscripts  denote  the  time 
steps,  while  solutions  for  the  first  and  the  second  time  steps  are  known. 

Equation  (27)  can  be  imposed  in  a  similar  fashion.  This  condition, 
within  the  limit  of  small  perturbation  assumption,  will  ensure  flow  just  above 
and  below  the  branch  cut  behaves  as  an  unsteady  vortex  sheet,  and  in  the 
steady  case  this  vortex  sheet  will  vanish  eventually. 

3.  NUMERICAL  RESULTS 

As  mentioned  earlier,  the  Galerkiii  approach  was  found  to  be  invalid 
for  computing  critical  transonic  flows  and  hence  abandoned.  As  an  alterna¬ 
tive,  theleast  squares  approach  was  pursued  and  applied  to  compute  both 
steady  and  unsteady  transonic  flows.  Results  obtained  thus  far  indicate  that 
the  least  squares  approach  is  generally  very  satisfactory.  Computations 
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performed  using  the  unsteady  code  are  presented  and  discussed  in  the  follow¬ 
ing. 

Steady  Transonic  flow  over  a  6  Thick  Circular  Arc 

In  order  to  debug  the  unsteady  code  and  study  the  feasibility  of  the  pre¬ 
sent  approach,  steady  state  solutions  for  flow  over  a  6%  thick  circular  arc 
airfoil  were  first  computed  and  compared  to  those  obtained  previously  by 
solving  the  steady  transonic  equation.  The  mesh  used  in  the  computations  is 
shown  in  Fig.  10,  where  the  typical  time  step  used  is  At  =  1. 

The  case  of  M  =  0.861  (barely  critical)  was  computed  first.  Conver- 

00 

gent  solutions  was  obtainable  whether  or  not  an  improper  downwind  influence 
un  the  solution  at  an  upwind  station  in  the  supersonic  pocket  is  discarded  or  not. 
The  two  solutions  are  unnoticeably  different.  This  is  obviously  because,  for 
=  0.861,  the  local  supersonic  region  is  very  small  and  therefore  the  effect 
of  the  improper  upwind  influence  is  insignificant.  Shown  in  Fig.  11  are  results 
eliminating  the  downwind  influence  and  compared  with  experimental  results 
(Ref.  13). 

For  cases  with  higher  freestream  Mach  number,  however,  numerical 
solutions  corresponding  to  with  and  without  downwind  influence  in  the  spatial 
direction  are  definitely  different.  Solution  with  improper  downwind  influence 
removed  produces  a  shock  wave,  although  somewhat  smeared;  while  the  other 
case  yields  a  shockless  solution  with  flow  expanded  and  recompressed  smoothly 
throughout  the  supersonic  pocket.  The  latter  results  had  also  been  observed 
in  the  steady  equation  formulation,  when  downwind  influence  was  retained  and 
Newton -Raphson’s  method  was  used  to  solve  the  nonlinear  algebraic  equations. 
Nevertheless,  solution  of  the  latter  type  (shockless  flow)  is  considered  to  be 
hydrodynamically  incorrect  due  to  considerations  of  entropy  production  through 
shocks.  For  this  reason,  the  scheme  with  improper  downwind  influence  re¬ 
moved  in  the  supersonic  region  is  adapted.  Results  for  the  cases  of  M  =  0.909 

oo 

and  =  0.966  are  also  shown  in  Fig.  11  and  compared  with  experimental  data 
due  to  Knechtel  (Ref.  13).  It  is  seen  that  the  shock  position  for  the  two  cases  is 


't 
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somewhat  behind  that  found  in  experiment,  a  fact  also  observed  in  finite 
difference  calculations.  In  our  numerical  scheme,  shock  continuity  has  not 
been  treated  as  such,  and  the  shock  was  smeared  at  least  within  the  width  of 
the  element  where  the  shock  is  supposed  to  be  located.  A  simple  way  to  obtain 
a  belter  resolution  is  to  refine  the  mesh  in  the  recompression  region. 

On  attempting  to  refine  the  solutions  for  M  =  0.909  and  0.966,  a  mesh 

00 

layout  depicted  in  Fig.  12  was  used.  The  results  for  the  case  of  M  =  0.909 

00 

agree  very  well  with  those  obtained  from  using  the  coarser  mesh,  except  for 

flow  after  the  shock,  where  recompression  is  stronger  for  the  finer  mesh. 

This  is  expected,  when  the  mesh  is  refined  the  (smeared)  shock  is  confined 

within  a  smaller  region,  and  better  resolution  should  result.  For  the  case 

of  M  =  0.966  the  shock  was  also  found  to  be  closer  to  the  trailing  edge,  com¬ 
oo 

pared  with  the  solution  obtained  using  the  coarser  mesh.  The  present  results 
appear  to  agree  with  Murman's  recent  finding  (Ref.  14)  about  embedded  shock 
over  airfoil.  He  observed  that  when  FCR  (fully  conservative  relaxation)  is  used, 
the  solution  gives  stronger  shocks  that  are  farther  aft  on  the  airfoil  compared 
with  that  using  the  NCR  (not  fully  conservative  relaxation).  For  reference, 
his  discussion  on  why  NCR  agrees  better  with  experimental  results  than  FCR 
may  be  consulted. 

Unsteady  Transonic  Flow  over  a  NACA  64  A006  Airfoil 

To  demonstrate  the  applicability  of  the  present  approach,  several  cases 
of  unsteady  transonic  flow  over  a  NACA  64  A006  airfoil  were  studied.  The 
airfoil  is  assumed  at  zero  incidence  but  with  a  quarter  chord  flap  executing 
harmonic  motion.  The  cases  studied  include 


M 

oo 

=  0.794, 

k  =  0.064 

M 

00 

=  0.804, 

k  =  0.253 

M 

00 

=  0.901, 

k  =  0.057 

M 

00 

=  0.903, 

k  =  0.228 
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The  element  mesh  used  is  shown  in  Fig.  13,  consisting  of  165  elements  with 

198  nodes.  For  each  cycle  of  oscillation,  16  time  steps  were  used  in  the 

computations.  For  the  two  supercritical  cases,  M  =  0.901  and  0.903,  unlike 

00 

their  steady  counterparts,  the  special  "one-sided  assembly"  scheme  was  not 
necessary  in  the  computations,  and  harmonic  solutions  are  apparently  obtained 
in  two  cycles  of  computations. 

Numerical  results  including  the  unsteady  pressure  jump  at  typical 
stations  (x/c  =  0.725  and  0.775)  and  the  amplitude  of  the  unsteady  pressure 
difference  over  the  airfoil  are  depicted  and  compared  in  Figs.  14  through  25. 

In  these  figures,  the  pressure  is  normalized  with  respect  to  the  angle  8,  the 
amplitude  of  flap  oscillation.  The  experimental  data  are  reported  by  Tijdeman 
and  Schippers  (Ref.  16).  In  general  the  agreement  for  predicted  results  and 
experimental  data  is  very  good.  In  particular,  results  for  the  typical  stations 
(x/c  =  0.725  and  0.775)  near  the  hinge  point  show  excellent  agreement  in  both 
magnitude  and  phase  angle.  The  amplitude  of  the  unsteady  pressure  difference 
over  the  airfoil  is  al;.u  depicted  with  the  comparisons  made.  Again,  the  experi¬ 
mental  data  have  been  converted  to  yield  the  value  of  amplitude.  The  dashed 
line  in  the  figure  corresponds  to  results  obtained  earlier  by  an  approach  de¬ 
veloped  for  NASA-Langley  (Ref .  6),  assuming  the  unsteadiness  to  be  a  small 
perturbation  to  the  mean  steady  flow.  Note  that  the  results  obtained  by 
the  two  theoretical  approaches  agree  very  well  over  most  parts  of  the  airfoil; 
and  agreement  between  present  results  and  experimental  data  is  very  good  for 
the  aft  half  of  the  airfoil.  Close  to  the  leading  edge,  both  finite  element  solu¬ 
tions  appear  to  under -predict  the  amplitude  of  unsteady  pressure  jump.  Many 
factors  could  attribute  to  this  fact.  The  small  perturbation  assumption  cur¬ 
rently  used  requires  that  flow  be  only  slightly  perturbed  from  uniform  flow. 
Consequently,  the  flow  behavior  in  the  vicinity  of  a  leading  edge  stagnation 
point  cannot  be  correctly  represented  with  such  an  assumption.  The  coarse¬ 
ness  of  the  present  mesh  may  also  cause  inaccuracy,  especially  near  the 
leading  edge.  Furthermore,  the  assumption  of  no  disturbance  along  the  outer 
boundary  of  the  mesh  as  used  in  the  present  computations  is  considered  to  have 
a  tendency  to  suppress  the  development  of  a  pressure  jump,  especially  for 
the  leading  edge  region,  where  no  oscillation  is  taking  place. 
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Figure  18  -  Time  History  of  Unsteady  Pressure  at  x/c  =  0.775 
(M  =  0.804,  k  =  0.253) 


It  is  seen  from  the  predicted  results  that  the  present  approach  is  suited 
for  computing  unsteady  transonic  flows  and  steady  flow  as  well.  For  an  airfoil 
undergoing  harmonic  motion,  oscillatory  solution  is  obtainable  in  two  to  three 
cycles,  and  requires  only  moderate  computation  time.  For  the  mesh  presently 
used,  marching  in  time  one  step  forward  requires  approximately  40  seconds 
CPU  time  on  a  Univac  1108  computer,  and  each  case  requires  about  10.5 
minutes  for  completing  computations  for  one  cycle,  with  16  time  steps.  As 
the  present  computer  program  has  not  been  optimized  yet,  some  reduction  in 
computation  time  is  still  possible. 


SECTION  IV 


CONCLUSIONS  AND  RECOMMENDATIONS 

The  work  presented  in  this  report  represents  an  initial  attempt  toward 
developing  a  general  numerical  method  for  computing  steady  and  unsteady 
transonic  flows  over  thin  airfoils  based  on  small  disturbance  theory.  Unlike 
other  existing  approaches  for  calculating  unsteady  transonic  flows,  the  pre¬ 
sent  numerical  algorithm  solves  directly  the  unsteady  transonic  equation 
with  both  time  derivative  terms  retained  in  the  computations.  Thus  the  pre¬ 
sent  algorithm  can  be  applied  to  compute  a  much  wider  class  of  transonic 
flow  problems,  including  steady,  oscillatory  or  transient  solutions,  either 
with  or  without  angle  of  attack.  For  oscillatory  flow,  no  assumption  is  made 
regarding  the  oscillating  frequencies,  nor  is  the  unsteady  perturbation  neces¬ 
sarily  to  be  small  compared  to  the  mean  steady  solution. 

The  present  numerical  algorithm  is  developed  using  the  finite  element 
concept  in  conjunction  with  the  least  squares  method  of  weighted  residuals 
applied  to  both  space  and  time.  Also,  schemes  for  handling  mixed  flows  and 
problems  involving  an  infinite  domain  have  been  developed  and  implemented 
in  the  present  numerical  method.  These  newly  developed  schemes,  together 
with  the  advantages  offered  by  the  finite  element  technique  which  includes 
great  flexibility  in  mesh  arrangement,  ease  with  implementing  higher  order 
approximations,  and  its  effectiveness  in  treating  complex  geometry  and 
boundary  conditions,  etc.,  make  the  present  approach  feasible  and  should 
prove  to  be  definite  assets  in  computing  more  complicated  transonic  flow 
problems . 

The  developed  procedures  have  been  applied  to  solve  several  problems 
of  steady  and  unsteady  transonic  flows  over  thin  airfoils,  and  results  obtained 
compare  generally  very  well  with  experimental  data  and  those  obtained  by 
other  numerical  techniques.  Regarding  computational  efficiency,  only  moderate 


computation  time  is  needed,  even  the  most  time  consuming  case  requires 
only  about  ten  minutes  CPU  time  on  a  CDC  6600  computer.  These  explora¬ 
tory  studies  well  indicate  the  feasibility  and  applicability  of  the  finite  ele¬ 
ment  approach  for  computing  steady  and  unsteady  transonic  flows.  However, 
as  the  effort  to  date  has  concentrated  on  investigating  the  best  suited  finite 
element  approach,  little  consideration  has  been  given  to  optimize  the  compu¬ 
tational  efficiency  and  to  certain  details,  such  as  more  accurate  treatment  of 
the  leading  edge  singularity  and  unsteady  effects  along  the  outer  boundary. 
These  refinements  are  apparently  desiraDle  and  certainly  could  be  accom¬ 
plished.  In  addition  to  these  improvements,  it  is  recommended  that  further 
study  be  pursued  in  the  following  areas: 

a.  Some  studies  need  to  be  done  on  the  effects  of  "non-conforming"  versus 
"conforming"  elements.  As  is  known,  the  elements  presently  used  are 

continuous  only  at  nodal  points,  but  not  across  element  boundaries, 
which  apparently  violates  (^continuity  between  elements  as  nor¬ 
mally  required.  Although  our  experiences  and  numerical  results 
indicate  that  the  use  of  "non -conforming  "  elements  do  not  cause 
problems  with  convergence,  as  concluded  in  plate  bending  analysis, 
some  studies  on  the  effects  of  "non -conforming "  versus  "conforming" 
elements  are  apparently  highly  desirable  and  useful,  considering  the 
fact  that  the  present  problem  is  nonlinear  and  involves  equations  of 
the  mixed  elliptic -hyperbol ic  type.  "Conforming"  elements  do  exist 
and  could  be  readily  implemented.  This  study  could  be  conducted 
concurrently  with  the  development  of  a  more  sophisticated  scheme  for 
handling  shock  waves,  such  as  shock-fitting. 

b.  The  procedures  developed  in  this  study  could  be  extended  without 
too  much  effort  to  treat  problems  involving  multi -element  airfoils 
and  the  cascade  problem.  As  techniques  for  handling  shock  waves 
(by  smearing)  and  Kutta  conditions  in  the  wake  have  already  been 
developed  and  appear  to  work  well  for  a  single  airfoil,  the  antici¬ 
pated  major  effort  in  solving  the  problem  of  multi -element  airfoils 
and  cascades  would  be  to  develop  or  implement  some  efficient 
scheme  for  solving  a  large  system  of  algebraic  equations,  such  as  an 
out-of-core  equation  solver.  Alternatively,  a  scheme  based  on  the 
concept  of  line  relaxation,  as  extensively  used  in  finite  difference 
techniques,  could  be  developed  to  cope  with  the  problem. 

c.  It  is  desirable  and  apparently  feasible  to  develop  a  finite  element 
algorithm  to  solve  the  nonlinear  full  potential  equation,  as  the  versa¬ 
tility  of  the  method  can  be  more  effectively  utilized  for  such  a  prob¬ 
lem.  Another  possibility  is  to  solve  the  two-dimensional  unsteady 
Euler  equations  directly  and  at  the  end  include  also  the  viscous  effects. 
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F  inite  difference  relaxation  techniques  for  solving  these  equations 
have  already  been  developed;  however,  tremendous  computation 
time  is  usually  required  even  for  solving  a  two-dimensional  prob¬ 
lem.  In  the  light  of  fully  utilizing  the  advantages  offered  by  the 
finite  element  method,  such  as  great  flexibility  in  mesh  arrange¬ 
ment,  ease  with  implementing  higher  order  approximations,  and 
its  effectiveness  in  treating  complex  geometry  and  boundary  con¬ 
ditions,  the  method  could  lead  to  higher  computational  efficiency 
than  other  existing  numerical  techniques.  Through  these  studies, 
the  feasibility  of  using  finite  element  methods  to  solve  the  three- 
dimensional  and  more  difficult  transonic  flow  problems  can  also 
be  assessed. 
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